###Replication data for "Does inequality beget inequality? Experimental tests of the prediction that inequality increases system justification motivation"
### Kris-Stella Trump and Ariel White
### December 2017
### This replication code for Additional Study 2 presented in Supplemental Information
### Code tested in R version 3.4.1

#Create log file
sink(file="Add_Study2_analysis_code_log.txt")

#Packages
install.packages("psych")

require(psych)

#Loading data 
#Don't forget to set working directory
mturk.data = read.csv("Congress MTurk data.csv")
online.data = read.csv("Congress Online sample data.csv")
names(mturk.data)
names(online.data)

#Indicator of recruitment source
mturk.data$mturk <- c(1)
online.data$mturk <- c(0)
#Combine data sources
cdata <- merge(mturk.data, online.data, all.x=T, all.y=T)


#### data prep
#For each experimental condition, remove the "NA" values and replace them with 0
cdata[is.na(cdata$cong.low)==T,"cong.low"] <- 0
cdata[is.na(cdata$cong.high)==T,"cong.high"] <- 0
cdata[is.na(cdata$inesc.low)==T,"inesc.low"] <- 0
cdata[is.na(cdata$inesc.high)==T,"inesc.high"] <- 0

#Create control condition indicator 
cdata$control <- ifelse(cdata$cong.low==0 & cdata$cong.high==0 & cdata$inesc.low==0 & cdata$inesc.high==0, 1,0)


#Code dependent variables

##SJT outcome 
class(cdata$sjt_1)
summary(cdata$sjt_1)
#Reverse code items 3 and 7
cdata$sjt_3r <- 10-cdata$sjt_3
cdata$sjt_7r <- 10-cdata$sjt_7
#Combined sjt score
cdata$sjt <- (cdata$sjt_1 + cdata$sjt_2 + cdata$sjt_3r + cdata$sjt_4 + cdata$sjt_5 + cdata$sjt_6 + cdata$sjt_7r + cdata$sjt_8)/8
summary(cdata$sjt)


##Gender outcome 
cdata$gender_orig <- cdata$gender #(keep the original MTurk "gender" column)
#Gender dv - recode so higher numbers indicate more justification of differences
cdata$gender.dv_1r <- 7-cdata$gender.dv_1
cdata$gender.dv_3r <- 7-cdata$gender.dv_3
cdata$gender.dv_4r <- 7-cdata$gender.dv_4
cdata$gender.dv_5r <- 7-cdata$gender.dv_5

cdata$gender <- (cdata$gender.dv_1r + cdata$gender.dv_2 + cdata$gender.dv_3r + cdata$gender.dv_4r + cdata$gender.dv_5r)/5
summary(cdata$gender)

##Trust DV
#All in one direction, recode so higher number is more trust
cdata$gss.conf <- (cdata$gss.conf.1_1 + cdata$gss.conf.1_2 + cdata$gss.conf.1_3 + cdata$gss.conf.1_4 + cdata$gss.conf.1_5 + cdata$gss.conf.1_6)/6

cdata$trust <- 4-cdata$gss.conf
summary(cdata$trust)

#### Break down cdata into the three experiments with different DV's
congress <- cdata[cdata$cong.low==1 | cdata$cong.high==1,]
inesc <- cdata[cdata$inesc.low==1 | cdata$inesc.high==1,]
control <- cdata[cdata$control==1,]

##where did participants find congress study
sum(congress$Q63==1, na.rm=T) #Harvard
sum(congress$Q63==2, na.rm=T) #SocialPsych.org
sum(congress$Q63==4, na.rm=T) #Social Psychology on the Net


#### Basic descriptive stats:
##all data
cdata$yearborn <- 1995-cdata$born
cdata$age <- 2014-cdata$yearborn
summary(cdata$age) 
table(cdata$gender_orig) 
summary(cdata$educ) 
length(unique(cdata$state))
summary(cdata$hhsize) 

##Mturk sample:
#names(mturk.data)
mturk.data$yearborn <- 1995-mturk.data$born
mturk.data$age <- 2014-mturk.data$yearborn
summary(mturk.data$age) 
table(mturk.data$gender) 
summary(mturk.data$educ) 
length(unique(mturk.data$state))
summary(mturk.data$hhsize) 

#online sample from elsewhere:
names(online.data)
online.data$yearborn <- 1995-online.data$born
online.data$age <- 2014-online.data$yearborn
summary(online.data$age) 
table(online.data$gender)
summary(online.data$educ) 
length(unique(online.data$state)) 
summary(online.data$hhsize)


#Internal consistency
#There are NA's in this data, but that is individuals who were in a different condition. Minimal NA's due to non-response.
require(psych)
#Congress study
alpha(congress[c("gss.conf.1_1", "gss.conf.1_2", "gss.conf.1_3", "gss.conf.1_4", "gss.conf.1_5", "gss.conf.1_6")])
alpha(congress[c("gender.dv_1r", "gender.dv_2", "gender.dv_3r", "gender.dv_4r","gender.dv_5r")])
alpha(congress[c("sjt_1", "sjt_2", "sjt_3r", "sjt_4", "sjt_5", "sjt_6", "sjt_7r","sjt_8")])

alpha(inesc[c("gss.conf.1_1", "gss.conf.1_2", "gss.conf.1_3", "gss.conf.1_4", "gss.conf.1_5", "gss.conf.1_6")])
alpha(inesc[c("gender.dv_1r", "gender.dv_2", "gender.dv_3r", "gender.dv_4r","gender.dv_5r")])
alpha(inesc[c("sjt_1", "sjt_2", "sjt_3r", "sjt_4", "sjt_5", "sjt_6", "sjt_7r","sjt_8")])




###Results analysis

###Compare control condition to the 'low ineq' and 'low inesc' conditions
t.test(congress[congress$cong.low==1,"sjt"],control[,"sjt"])
t.test(inesc[inesc$inesc.low==1,"sjt"],control[,"sjt"])
t.test(congress[congress$cong.low==1,"gender"],control[,"gender"])
t.test(inesc[inesc$inesc.low==1,"gender"],control[,"gender"])
t.test(congress[congress$cong.low==1,"trust"],control[,"trust"])
t.test(inesc[inesc$inesc.low==1,"trust"],control[,"trust"])
###Compare control condition to the 'high ineq' and 'high inesc' conditions
t.test(congress[congress$cong.high==1,"sjt"],control[,"sjt"])
t.test(inesc[inesc$inesc.high==1,"sjt"],control[,"sjt"])
t.test(congress[congress$cong.high==1,"gender"],control[,"gender"])
t.test(inesc[inesc$inesc.high==1,"gender"],control[,"gender"])
t.test(congress[congress$cong.high==1,"trust"],control[,"trust"])
t.test(inesc[inesc$inesc.high==1,"trust"],control[,"trust"])


### Congress treatment
t.test(congress[congress$cong.low==1,"sjt"],congress[congress$cong.high==1,"sjt"])
t.test(congress[congress$cong.low==1,"gender"],congress[congress$cong.high==1,"gender"])
t.test(congress[congress$cong.low==1,"trust"],congress[congress$cong.high==1,"trust"]) 


#Check whether source of participant made a difference
sjtreg <- lm(formula=sjt ~ cong.low*mturk, data=congress)
summary(sjtreg)
genderreg <- lm(formula=gender ~ cong.low*mturk, data=congress)
summary(genderreg)
trustreg <- lm(formula=trust ~ cong.low*mturk, data=congress)
summary(trustreg)

#Check whether trust in congress in particular is impacted
t.test(congress[congress$cong.low==1,"gss.conf.1_6"],congress[congress$cong.high==1,"gss.conf.1_6"])

### Inescapability treatment
t.test(inesc[inesc$inesc.low==1,"sjt"],inesc[inesc$inesc.high==1,"sjt"])
t.test(inesc[inesc$inesc.low==1,"gender"],inesc[inesc$inesc.high==1,"gender"])
t.test(inesc[inesc$inesc.low==1,"trust"],inesc[inesc$inesc.high==1,"trust"])

#put everything on a comparable scale (compress to 0-1)
summary(congress$sjt/9); summary(congress$gender/6); summary(congress$trust/4)
#SJT is 1-9, gender 1-6, trust 1-4
congress$sjtnew <- congress$sjt/9
congress$gendernew <- congress$gender/6
congress$trustnew <- congress$trust/4
inesc$sjtnew <- inesc$sjt/9
inesc$gendernew <- inesc$gender/6
inesc$trustnew <- inesc$trust/4

### Inescapability treatment on standardized scale
t.test(inesc[inesc$inesc.low==1,"sjtnew"],inesc[inesc$inesc.high==1,"sjtnew"])
t.test(inesc[inesc$inesc.low==1,"gendernew"],inesc[inesc$inesc.high==1,"gendernew"])
t.test(inesc[inesc$inesc.low==1,"trustnew"],inesc[inesc$inesc.high==1,"trustnew"])

### Congress treatment on standardized scale
t.test(congress[congress$cong.low==1,"sjtnew"],congress[congress$cong.high==1,"sjtnew"])
t.test(congress[congress$cong.low==1,"gendernew"],congress[congress$cong.high==1,"gendernew"])
t.test(congress[congress$cong.low==1,"trustnew"],congress[congress$cong.high==1,"trustnew"]) 


### time spent on treatments - congress
timec <- c(cdata$cong.low.t_3, cdata$cong.hight_3)
summary(timec) #median 19 mean 39
t.test(cdata$cong.low.t_3, cdata$cong.hight_3) #no difference in time spent by condition
sum(timec>5, na.rm=T) #253 of 266 people spend more than 5 seconds

### time spent on treatments - inesc
timec <- c(inesc$inesc.l.ti_3, inesc$inesc.h.ti_3)
summary(timec) #median 27 mean 32
t.test(inesc$inesc.l.ti_3, inesc$inesc.h.ti_3) #no difference in time spent by condition
sum(timec>5, na.rm=T) #238 of 240 people spend more than 5 seconds

###Figure A6 in Supplemental Information 

##congress/various outcomes
cong_grid <- as.data.frame(cbind(rep(NA,3), rep(NA,3), rep(NA,3), rep(NA,3)))
cong_grid$V1 <- c("Gender \nJustification", "Institutional \nTrust", "System \nJustification")
colnames(cong_grid)<- c("outcome", "diff.est", "conf.low", "conf.high")

cong_gen<- t.test(congress[congress$cong.high==1,"gendernew"],congress[congress$cong.low==1,"gendernew"])
cong_gen$conf.int[2]
cong_grid[1, 2] <- cong_gen$estimate[1]- cong_gen$estimate[2]
cong_grid[1, 3] <- cong_gen$conf.int[1]
cong_grid[1, 4] <- cong_gen$conf.int[2]

cong_trust <- t.test(congress[congress$cong.high==1,"trustnew"],congress[congress$cong.low==1,"trustnew"])
cong_grid[2, 2] <- cong_trust$estimate[1]- cong_trust$estimate[2]
cong_grid[2, 3] <- cong_trust$conf.int[1]
cong_grid[2, 4] <- cong_trust$conf.int[2]
cong_grid$ruler <- c(1,2,3)

cong_sjt <- t.test(congress[congress$cong.high==1,"sjtnew",],congress[congress$cong.low==1,"sjtnew"])
cong_sjt$conf.int[1]
cong_grid[3, 2] <- cong_sjt$estimate[1]- cong_sjt$estimate[2]
cong_grid[3, 3] <- cong_sjt$conf.int[1]
cong_grid[3, 4] <- cong_sjt$conf.int[2]

pdf("congress_differenceests.pdf")
par(las=1, oma=c(1, 4, 1, 0))
plot(NA,NA, main="Effect of \nCongress Treatment", ylim=c(.5, 3.5), xlim=c(-0.12,0.12), yaxt='n', 
     xlab= "", ylab="", cex.main=1.3)
points(cong_grid$diff.est, cong_grid$ruler, pch=18, cex=2)
mtext("Estimated Difference Between \nHigh and Low Inequality Conditions", cex=1.3, side=1, line=4)
for (i in 1:3){
  segments(cong_grid$conf.low[i], cong_grid$ruler[i], cong_grid$conf.high[i], cong_grid$ruler[i])}
axis(side=2, at=cong_grid$ruler, labels = cong_grid$outcome, cex.axis=1.3)
abline(v=0, lty=3, col="gray")
dev.off()

